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The interaction between two bright solitons in a dipolar Bose-Einstein condensate (BEC) has 
been investigated aiming at finding the regimes when they form a stable bound state, known as 
soliton molecule. To study soliton interactions in BEC we employed a method similar to that used 
in experimental investigation of the interaction between solitons in optical fibers. The idea consists 
in creating two solitons at some spatial separation from each other at initial time to, and then 
measuring the distance between them at a later time ti > to- Depending on whether the distance 
between solitons has increased, decreased or remained unchanged, compared to its initial value at to, 
we conclude that soliton interaction was repulsive, attractive or neutral, respectively. We propose 
an experimentally viable method for estimating the binding energy of a soliton molecule, based 
on its dissociation at critical soliton velocity. Our theoretical analysis is based on the variational 
approach, which appears to be quite accurate in describing the properties of soliton molecules in 
dipolar BEC, as reflected in good agreement between the analytical and numerical results. 


I. INTRODUCTION 


The interaction of solitons has been a subject of great interest right from the beginning of their early investi¬ 
gations pj. New fundamental features of soliton interactions are still being discovered, in fact the existence of a 
phase-dependent spatial jump in the trajectories of two colliding matter-wave solitons, reported in a recent exper¬ 
iment Q, is just one example to be mentioned. Apart from its scientific importance soliton interactions have also 
a practical importance. For instance, interaction between optical solitons sets the limit on the rate of information 
transfer in fiber optic communication systems Q. Due to important applications, soliton interactions are extensively 
studied, both theoretically and experimentally, in optical fibers SH3, photonic crystals Q and plasmas Recent 
experimental studies have shown that in addition to interactions between neighboring optical solitons in close prox¬ 
imity, there exists also a long range interaction between them 0- The generation of spatially separated coherent 
matter-wave packets and_ their subsequent interaction constitute the basic phenomena in the operation of modern 
atomic interferometers El working in the solitonic regime fl2l . IT3ll where the fringe visibility is significantly increased 
compared with an ordinary atomic cloud, as demonstrated in |12j. 

Solitons have been experimentally observed in many areas of physics, including the Bose-Einstein condensates 
(BEC) [110 . The experiments with solitons in BEC reported so far were concerned with the creation of solitons 
and studying their collective dynamics. Regarding the type of interaction between matter-wave solitons a conjecture 
was made from the behavior of neighboring solitons in a soliton train El- Meanwhile, it would be interesting to 
explore systematically the interaction between two matter-wave solitons with varying spatial separation and relative 
phase. Recent progress in controlled creation and manipulation of matter-wave solitons in BEC [11110 indicates 
that such experiments on soliton interactions are now within the scope of current technology. A key role belongs 
to a minimally destructive polarization phase-contrast imaging technique 0, that allows to make multiple imaging 
of the soliton pair during a single experimental run, as reported recently with regard to phase-dependent collision 
of two matter-wave solitons Q. An essentially new method reported in Ref. [16] for controlled (i.e. deterministic 
in both soliton position and momentum) creation of matter-wave bright solitons and soliton pairs without the use 
of Feschbach resonances opens new perspectives for investigation of soliton interactions in BEC with unprecedented 
accuracy. 

Experimental realization of chromium BEC with long range dipole-dipole atomic interactions [19 ) has opened new 
direction in the physics of ultra-cold quantum gases. Subsequently two other species with strong dipolar interactions, 
namely dysprosium [2Cj and erbium 21 [, were Bose-condensed. The principal difference of chromium condensates from 
the alkali atom condensates is that, 52 Cr has a large permanent magnetic dipole moment d = 6 / ib , where /is = eh^rrie 
is the Bohr magneton. Since the dipole-dipole force is proportional to the square of the magnetic moment, the dipolar 
interactions in chromium condensate is a factor of 36 times stronger than in alkali atom condensates, like 8 'Rb 
(d = 1 fis)- Similar arguments pertain also for other dipolar quantum gases, 164 Dy ( d = 10 hb) and 168 Er ( d = 7 hb)- 

In this work we study, by means of variational approximation (VA) and numerical simulations, the interaction 
between two bright solitons in a dipolar BEC. We employ a strategy similar to that used in the experimental inves¬ 
tigation of the interaction forces between fiber optic solitons [HI. |22|. Following that idea in numerical experiments 
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we create two bright solitons at some initial spatial separation, then give the pair a chance to evolve for some period 
of time, and finally measure the distance between the solitons when the evolution time has elapsed. Depending on 
whether the distance between the solitons has increased, decreased or remained unchanged, compared to its initial 
value, we conclude about the type of soliton interaction as being repulsive, attractive or neutral, respectively. 

There is a qualitative difference between solitons in dipolar and non-dipolar media. Specifically, two anti-phase 
solitons in dipolar media attract each-other at large separation and repel at short separation. Due to this property 
they can form stable bound states with non-zero binding energy, whereas in non-dipolar media they always repel 
and never form stable bound state. The possibility of molecular type of interaction between solitons in dipolar BEC 
moving in two neighboring wave- guid es was shown in [23}. The existence of stable multi-soliton structures in 2D 
dipolar BEC was also reported in [24j . 

Our main objective in this work is to find the conditions when two interacting solitons in the same quasi-ID 
waveguide form a stable bound state, which can be considered as a basic matter-wave soliton molecule. When 
a stable bound state of two solitons has been realized, we characterize the soliton molecule by its bond length 
and binding energy. Our work distinguishes itself from other relevant publications in that, we use the VA with a 
Gauss-Hermite ansatz and analytically tractable function of non-locality (response function), which allows to describe 
essential features of soliton molecules in dipolar BEC. Moreover, we provide detailed comparison of predictions of VA 
with the results of numerical simulations of the Gross-Pitaevskii equation. 

The paper is organized as follows. In the next section II we introduce the governing equation and develop the 
VA for the dynamics of soliton molecules in dipolar BEC. In Sec. Ill we use an optimization procedure to find the 
shape of a soliton molecule and validate the VA by comparing the analytical predictions with the results of numerical 
simulations. In Sec. IV we present such an important characteristic of a soliton molecule as its binding energy. In 
Sec. V we reveal the character of soliton interactions in dipolar BEC using the method borrowed from the field of 
fiber optic solitons. In Sec. VI we summarize our findings. 


II. THE GOVERNING EQUATION AND VARIATIONAL APPROACH 

From the viewpoint of theoretical description, matter-wave solitons in BEC and optical solitons in fibers are similar. 
The mean field Gross-Pitaevskii equation (GPE) for the dynamics of BEC and the nonlinear Schrodinger equation 
for propagation of optical solitons in fibers have formal analogy. Similarity of the basic equations has been fruitful in 
transferring many ideas from nonlinear optics to the field of matter-waves [25} . In this paper we transfer one more 
idea, concerning soliton interactions, from the field of fiber optics into the field of BEC. 

We shall consider the one dimensional GPE by taking into account both the local and nonlocal nonlinearities, which 
account for the usual contact interactions between atoms, and the long-range dipole-dipole interactions 12(1 127} 
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where q = a s /|a s o| is the coefficient of contact interactions, controlled by the atomic s-wave scattering length a s , with 
a s o being its background value, g = ad/|a s o| is the coefficient of nonlinearity, responsible for the long - range dipolar 
atomic interactions, expressed via characteristic dipole length a <2 = god 2 m/(\2'i:‘h 2 ), with m,d being the mass and 
magnetic dipole moment of atoms, oriented along the x axis, go is the permeability of vacuum. Time and space are 
expressed in units of t 0 = wj) 1 and l 0 = \/h/(mui±_), respectively, with w_l being the frequency of radial confinement. 
The wave function is re-scaled as i/j = y / 2|a s o|'I' and normalized to the reduced number of atoms in the condensate 
N = \ijj(x)\ 2 dx, which is a conserved quantity of Eq. dT|). The following two models for the kernel (nonlocal 
response functions) are relevant to dipolar condensates confined to quasi-ID traps 

Ri(x) = (1 + 2x 2 ) exp(a; 2 ) erfc(|x|) — 2t:~ 1 ^ 2 \x\, (2) 

R 2 {x) = S 3 (x 2 + S 2 )~ 3/2 . (3) 

The former kernel was derived for the dipolar BEC using the single mode approximation [28} . while the latter 
containing a cutoff parameter d, was proposed in Ref. [26] and is more convenient for analytical treatment. Making 
use of the matching conditions 
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which requires 6 = 7 r -1 / 2 , one can take advantage of the simplicity of R, 2 (x) for the application of VA. The meaning 
of S is the effective size of the dipole. Actually, it takes value of the order of the transverse confinement length, which 
makes the model one-dimensional, and sets the unit length in Eq. ©■ Therefore, the choice of 5 = tt 1 / 2 » 0.56 is 
quite reasonable. In the limit x 5, where dipole-dipole interaction effects dominate the contact interaction effects, 
both response functions behave as ~ 1/x 3 . This justifies application of the kernel function R, 2 {x) for analytical 
treatment of dipolar effects in BEC. By comparing the graphics of these two response functions one can be convinced, 
that indeed I?i ( x ) and R 2 (x) match very closely [2(|. 

The Lagrangian density generating the Eq. © is 

00 

£ = - V>*tM + \\^\ 2 - |M 4 - !|V>(M)I 2 J R{x-t)\i>(t,t)\ 2 dt. (5) 


To study soliton interactions in dipolar BEC we need to develop the VA for a two-soliton molecule. To this end 
we employ a Gauss-Hermite (GH) trial function, which was successful in the description of soliton molecules in 
dispersion-managed optical fibers 1291 
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where A(t), a(t), b{t ) and (f>(t) are variational parameters, associated with the amplitude, width, chirp and phase, 
respectively. The norm N = J \ip(x)\ 2 dx = A 2 o 3 \/ : k/2 is proportional to the number of atoms in the condensate. For 
specified values of A and a, the waveform © can be well approximated by two anti-phase Gaussian functions with 
the amplitude Aq, width ao and half separation xq 
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2Aa 2 /tt 
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Substitution of the ansatz © and response function © into the Lagrangian density © and subsequent integration 
over the space variable x yields the averaged Lagrangian 
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is the confluent hypergeometric function [30], and z = 5 2 /(2a 2 ). 

The VA equation for the parameter a of the two-soliton molecule, which is proportional to the separation between 
solitons, can be derived from the Euler-Lagrange equations d/dt(dL/dvt) — dL/du = 0 for variational parameters 
v —> a, b, (j), using the averaged Lagrangian © 


a t t 


1 qN 
a 3 4\/27r a 2 


9$N 
4a/ 2 a 2 





( 10 ) 


The corresponding effective potential U(a) for the width is depicted in Fig. [T] (left panel). The analytic form of the 
potential 17(a), which can be found by integrating the right hand side of Eq. (fTOl) . is rather complicated and we do 
not show it here explicitly. The fixed point att = —dU(a)/da = 0 of this equation a 0 is associated with the stationary 
separation between center-of-mass positions of two solitons constituting the molecule Ao = 2xq = 4ao j\fn. At larger 
separation ( a > ao) the solitons attract each other (dU/da > 0), and at smaller separation (a < ao) they repel 
( dU/da < 0), therefore the effective potential U(a) has a property of molecular type. Fig. [T| (right panel) illustrates 
the shape of a two-soliton molecule, as predicted by VA, by two anti-phase Gaussian functions with parameters given 
in Eq. ©, and by optimization procedure, applied to GPE ©, described in the next subsection. 

When solitons of the molecule are placed at their equilibrium positions, they stay motionless, as shown Fig. [2] 
If solitons are slightly displaced and released, they perform small amplitude oscillations around their stationary 
separation. The dynamics of the molecule strongly depends on the initial phase difference between solitons. In 
particular, even slight deviation from anti-phase configuration leads to periodic exchange of atoms between solitons. 
At larger deviation the soliton molecule does not form. 
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FIG. 1: (Color online) Left panel: A molecular type potential associated with VA Eq. (1101) for different strengths of the contact 
interaction. Right panel: The shape of a two-soliton molecule in a pure dipolar BEC with N = 2, q = 0, g = 20 as predicted 
by VA with trial function © for A = 3.538, a = 0.553, N = 1.876 (red solid line), by two anti-phase Gaussian functions with 
parameters given in Eq. © (blue dashed line), and found from numerical optimization procedure, applied to GPE © (brown 
dash-dot line). The minimum of the effective potential U(a) is attained at a o = 0.53, and the equilibrium distance between 
center-of-mass positions of pulses predicted by Eq. © is Aq = 4ao/ \Je — 1-2, while the GPE optimization gives Aq — 1.3. 
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FIG. 2: (Color online) (a) Stable propagation of a two-soliton molecule composed of two anti-phase (ifi = n) Gaussian pulses 
with parameters Aq = 1.196, a o = 0.394, placed at their equilibrium positions xo = ±0.657. (b) When the solitons are slightly 
(by 20 %) displaced from equilibrium positions, they perform oscillations. Density plot \ij)\ 2 is obtained by numerical solution of 
the GPE ©. Dashed line corresponds to calculations according to VA equation (1101) . which shows increasing phase shift with 
respect to GPE. (c) Periodic exchange of atoms between two solitons when the initial phase difference is slightly decreased. 


The frequency of soliton oscillations near equilibrium state can be estimated from linearized version of Eq. m 
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For the period of oscillations near the stationary separation we have Tva = 2n/£lo ~ 1.54. The prediction of VA is in 
qualitative agreement with the result of numerical simulation of the GPE, Tqpe — 2.2 (see Fig. ©. In general, the 
VA provides fairly good description of the static and dynamic properties of the soliton molecule, while its waveform 
remains close to the selected trial function (©). The agreement between VA and GPE deteriorates at large separation 
between solitons, close to the dissociation point, when the trial function cannot be well approximated by two anti-phase 
Gaussian functions. 
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III. IMPROVING THE SHAPE OF THE SOLITON MOLECULE BY OPTIMIZATION PROCEDURE 

The VA provides approximate waveform of a soliton molecule. When a trial function with parameters, defined by 
stationary solution of VA equation, is assigned as initial condition for the GPE, small amplitude oscillations of the 
molecule’s shape and separation between pulses is observed. This implies that soliton molecule is in its excited state. 

For some precise parameter calculations, such as the binding energy of soliton molecules, a truly ground state 
should be employed. In Ref. [3l| an optimization strategy to find the stationary shape of a soliton molecule in 
dispersion-managed optical fibers was proposed. Below we extend this approach for soliton molecules in dipolar BEC. 
It is based on the Nelder-Mead (NM) nonlinear optimization procedure 1321 . which seeks to minimize an objective (or 
cost) function 


/ = 


An 


(IV’(uO)l - \ip{x,tx)\) 2 dx, 
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\ip{x, 0)| 2 dx, 


( 12 ) 
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is the initial waveform, composed of two anti-phase Gaussian functions, separated by a distance 2xo, and ip(x,ti) is 
the result of evolution of ip(x 1 0) for some period of time t = ti, according to GPE. Normalization factor Nq in Eq. 
m is introduced to avoid trivial solutions, in particular corresponding to Xq = 0. In general case minimization of 
the objective function can be performed with respect to variables ao and lo, since the amplitude Aq is fixed by the 
norm of the Gaussian. However, numerical experiments show that, VA predicted values of ao and Ao for a single 
soliton are quite accurate, and minimization only with respect to pulse separation Xq can produce stationary state 
of the molecule. The evolution time t\ can be estimated as half period of oscillation for the molecule t\ = 7 r/flo. 
Although the NM optimization procedure finds the minimum of the objective function (1121) for Gaussian functions 
with broad range of parameters, the convergence rate can be improved by selecting the initial waveform close to 
the stationary state. The VA can provide the waveform which is close to the stationary state. We find stationary 
pulse separation Xq and norm of the soliton molecule A from NM optimization procedure. The obtained results were 
confirmed by alternative method of Luus-Jaakola BSl. Our preference of these optimization methods is motivated 
by several advantages, such as simplicity of programming (since calculation of function derivatives is not required), 
high convergence rate and reliability and effectiveness in locating the global minimum of the objective function. 


IV. INTERACTION POTENTIAL AND BINDING ENERGY OF SOLITON MOLECULES 

The binding energy of a soliton molecule Eb can be defined as the amount of energy, which is required for the 
dissociation of the molecule into two separate individual free solitons, far away from each other. In numerical simula¬ 
tions using the GPE, the process of dissociation can be implemented by assigning an initial velocity to each soliton in 
opposite directions ip = ipie lvx + ip2e~ lvx . If the velocity is smaller than some critical value v < v cr , solitons perform 
oscillations around their stationary positions, otherwise the molecule disintegrates into individual solitons, travelling 
in opposite directions, as illustrated in Fig. [3] 

In a “particle in potential well” picture this situation corresponds to the escape of the particle from the potential 
well at critical kinetic energy. The critical velocity determines the binding energy of the molecule Eb ~ u 2 r /2. Fig. []] 
illustrates the potential of interactions between the two solitons, normalized to binding energy, as a function of 
distance between solitons in units of stationary separation xo■ To construct the potential U{x) we assign velocity to 
solitons and determine the maximal and minimal value of the separation, which correspond to right and left classical 
turning points of the oscillating particle in the potential well. Repeating these calculations for velocities in the range 
v G [0, v cr \ we construct the potential, shown in the left panel of Fig. 0] As expected, the bigger norm A (or number 
of atoms) of the molecule corresponds to stronger potential, connecting solitons. 

The critical velocity v cr , at which the molecule disintegrates into far separated individual solitons, is determined 
from GPE simulations, by setting in motion the two bound solitons in opposite directions, as shown in Fig. [I] In 
the experiment pushing the solitons in opposite directions can be realized by means of a laser beam, directed into the 
center of the molecule, as it was used to split the condensate in two halves Q. The intensity of the laser beam can 
be varied to give desired initial velocity to solitons. 

About the repulsive interaction between two anti-phase solitons the following comment is appropriate. As experi¬ 
mentally demonstrated in [2] and theoretically shown in [35j , two colliding wave packets exchange not only velocities 
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FIG. 3: (Color online) Dynamics of center-of-masses of two solitons, forming the molecule, when the two solitons are set in 
motion with different velocities in opposite directions. When the velocity is less than critical v < v cr = 0.91, solitons perform 
oscillations near equilibrium positions. At critical velocity the molecule dissociates into freely moving individual solitons (red 
solid lines). 




FIG. 4: (Color online) Left panel: The potential of interaction between two solitons, retrieved from numerical GPE simulation, 
for two values of the molecule’s norm. Similarity with the VA predicted potential U(a) in Fig. [l] is evident. Right panel: The 
stationary half separation between solitons of the molecule xo (blue dashed line) and its binding energy Eh (red solid line) as a 
function of the molecule’s norm N. Symbols correspond to values found from numerical simulations of the GPE Q. and lines 
are interpolating curves for visual convenience. 


(as classical particles do), but also their entire wave-functions (as quantum mechanical particles do via tunnel phe¬ 
nomenon). In our case of equal masses of two colliding solitons, the classical and quantum descriptions lead to the 
same result. Physically, the repulsive interaction of anti-phase matter wave solitons can be regarded as exchange of 
velocities of two colliding classical particles interacting via hard core potential. 


V. NUMERICAL SIMULATION OF THE TWO-SOLITON INTERACTIONS 

In order to study the character of interaction between two matter-wave solitons in numerical experiments we employ 
the idea similar to that used for optical solitons in fibers [Bj]. Initially at t = 0, two solitons either in-phase (f> = 0 or 
out-of-phase </> = 7r, are created at some distance Ao from each-other. At later time t = t\ > 0 the distance between 
solitons is measured again. If the solitons did not interact, the distance between them should not change with respect 
to its initial value Ai = Ao. If the interaction was attractive, the final measured distance should be less than the 
initial distance Ai < Ao. Finally if the interaction was repulsive, the final distance should be greater than the initial 
distance A l > Aq. The numerical experiment consists in repeating the above procedure for different values of the 
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initial distance Ao, starting from sufficiently large separation, greatly exceeding the width of the soliton then reaching 
small distances when solitons start to overlap. 

The result is presented in Fig. [5] as a plot of Ai (final separation) vs. Ao (initial separation). When the solitons, 
comprising the molecule, are out of phase (<j> = 7r) we see that at large separations the solitons do not interact 
(Ai ~ Ao, red curve is close to median), while at small distances they attract each other (Ai < Ao, red curve is 
below the median) until they reach a minimum (stationary) separation, where again Ai = Ao (red dot on the median). 
When solitons are placed at even smaller distances, they repel each other (Ai > Ao, red curve is above the median). 
That is, the two-soliton molecule behaves like a diatomic molecule. For in-phase solitons (<j> = 0) we see that solitons 
attract each other until their separation becomes comparable to the width of the soliton, then merge forming a wave 
packet, whose shape strongly oscillates. 



FIG. 5: (Color online) Left panel: The character of soliton interactions for anti-phase and in-phase solitons. Anti-phase solitons 
attract each-other at large separations and repel at small separations (red solid line). There exists a stationary separation, 
shown by a red dot, where attraction changes to repulsion. In-phase solitons always attract and collide (blue dashed line). 
Middle panel: Anti-phase solitons can form a stable soliton molecule at appropriate initial separation Ao — 1.3. Right panel: 
In-phase solitons always collide and do not form stable bound state. 


The numerical simulations are performed using realistic values of atom numbers and interaction parameters in 164 Dy 
for which m = 2.7 x 10 -25 kg, d = 10/is = 9.27 x 10 - 23 A m 2 , ad = ^d 2 m/ ( 127 r?i 2 ) ~ 7x 10 -9 m. The frequency 
of radial confinement w_l = 27r x 62 Hz provides radial oscillator length Iq — 1 A™ and unit of time to = 2.6 ms. For 
parameter values go = 20 and N = 2 used in numerical simulations we obtain the number of atoms in a two-soliton 
molecule Af = goNlo/{2aa) — 3000. The total number of atoms in the 164 Dy condensate was AT = 15000 (2Cj. 


VI. CONCLUSIONS 

We have studied the interaction between two bright solitons in a dipolar BEC and found the conditions when 
they form a stable bound state. It was revealed by numerical simulations of the governing nonlocal GPE and 
corresponding variational analysis, that two anti-phase solitons in dipolar condensates behave similarly to a diatomic 
molecule. Namely, they attract each-other at large separation and repel each-other at small separation. There exists 
a particular distance at which the two solitons remain motionless in their stationary state. Solitons in a weakly 
perturbed molecule perform small amplitude oscillations near the equilibrium position, the frequency of which is 
predicted quite accurately by the developed model. Two in-phase solitons when placed close to each-other always 
collide and do not form the bound state. The obtained results can be useful in further studies of the properties of 
multi-soliton bound states in dipolar BEC. 
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